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Using N-body simulations, we measure the power spectrum of the effective dark matter density field, which 
is defined through the modified Poisson equation in f{R) cosmologies. We find that when compared to the 
conventional dark matter power spectrum, the effective power spectrum deviates more significantly from the 
ACDM model. For models with fm = —10“^, the deviation can exceed 150% while the deviation of the 
conventional matter power spectrum is less than 50%. Even for models with fm = —10“®, for which the 
conventional matter power spectrum is very close to the ACDM prediction, the effective power spectrum shows 
sizeable deviations. Our results indicate that traditional analyses based on the dark matter density field may 
seriously underestimate the impact of f{R) gravity on galaxy clustering. We therefore suggest the use of the 
effective density field in such studies. In addition, based on our findings, we also discuss several possible 
methods of making use of the differences between the conventional and effective dark matter power spectra in 
f[R) gravity to discriminate the theory from the ACDM model. 
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I. INTRODUCTION 

It is well established that the Universe is currently under¬ 
going a period of accelerated expansion Ul. The standard 
paradigm holds that this accelerated expansion is caused by 
a non-zero cosmological constant. A much discussed alterna¬ 
tive is that a modification to general relativity could also ac¬ 
count for the observed acceleration. However, the dynamics 
of the solar system heavily constrain the form that any such 
modification may take. Any viable theory of gravity must 
be indistinguishable from standard general relativity in high 
density environments, such as the Solar System and the cen¬ 
ters of galaxies. A viable theory must also reproduce the ob¬ 
served expansion history, which is very well fit by the standard 
ACDM cosmology 1321 . 

Chameleon f{R) gravity is a class of models where the ef¬ 
fective potential of the scalar field is dependent upon the lo¬ 
cal environment 0. A fifth force is involved which can be 
"screened" by the local density field in very deep potential 
wells. In regions where the matter density is high and the po¬ 
tential well is very deep, the fifth force is usually screened 
and gravity behaves just like standard GR. However, in low 
density regions of space there is no such screening and the 
strength of gravity is enhanced by the fifth force. 

There are several observational effects due to such a modi¬ 
fication. However, in order to make competitive forecasts for 
constraining /(i?) gravity with current and future cosmolog¬ 
ical surveys, it will be necessary to study the clustering of 
galaxies and to produce mock galaxy catalogues from simula¬ 
tions in f{R) gravity. 

However, the clustering of galaxies in f{R) gravity is very 
complicated. The gravitational force produced by the dark 
matter field is mediated by the fifth force, which is no longer 
the same as that in standard gravity. The relationship between 
dark matter halos and the clustering of galaxies is not the same 
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in f{R) gravity compared to the standard model. Although 
the clustering of dark matter halos in f{R) gravity has al¬ 
ready been studied in the literature 161, it might be a risk to 
use the standard dark matter halos in semi-analytical galaxy 
formation models fTl to analyze the formation and clustering 
of galaxies. However, if we define an effective density field, 
the gravitational force produced by the effective density field 
could still have the same form as that in standard gravity. The 
clustering of galaxies in the effective dark matter density field 
could still follow that in the standard gravity. It is therefore 
more interesting to study the effective density field than the 
standard dark matter density field when analyzing the forma¬ 
tion and clustering of galaxies. 


In Ref Is) using the effective density field, we studied the 
properties of effective halos. We found that the relationships 
between the effective mass of a halo and its dynamical prop¬ 
erties closely resemble those in the ACDM cosmology. This 
is an interesting result. The aim of this paper is to further ex¬ 
tend this idea. We shall not only focus on halos but also on 
the entire density field. We shall study the power spectrum of 
the effective density field since in f{R) gravity it is closely 
related to the galaxy power spectrum which is one of the most 
important statistical measures of the clustering of galaxies in 
cosmology. 


This paper is organized as follows; In Sec.lIIjwe introduce 
the f{R) model studied in this work. In Sec. |III[ we discuss 
the linear and non-linear perturbation equations in f{R) grav¬ 
ity. We also present the technical details of our simulations. 
In Sec. IV we describe our method of measuring the power 
spectra of scalar fields and we also present the numerical re¬ 
sults on the effective power spectra. In Sec. |V] we summarize 
and conclude this work. 




2 


II. F(R) MODEL 


The /(i?) gravity model is defined with the four¬ 
dimensional modified Einstein-Hilbert action 

S=^ J d'^xV^[R + fiR)] + J , ( 1 ) 


where — SttG with G being Newton’s constant, g is the 
determinant of the metric is the Lagrangian den¬ 

sity for matter, and f{R) is an arbitrary algebraic function of 
the Ricci scalar curvature R 191 4201 (see Refs. 12111221 for re¬ 
views). 

In this work, we consider an f{R) model that can exactly 
reproduce the ACDM background expansion history l23l . 
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The indices in the above expression are given by 
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2 F 1 [a, b; c; z] is the hypergeometric function. When c > b > 
0, the hypergeometric function has the integral representation 


= r(6)P(c- i.) I 

(3) 

where r(a;) is the Euler gamma function. 2 A'i[a, b; c; z] is a 
real function that is well defined in the range —00 < z < 1. 
Hq is the present Hubble constant. is the matter density 
today and = 1 ~ R additional parameter that 
characterises the f{R) model. Eor the instability issue as dis¬ 
cussed in Ref. l24l . D should be constrained as H < 0. The 
model predicts a lower bound for the scalar curvature R across 
the Universe 


R G (4A, + 00 ) , (4) 

where 

A = • (5) 


III. N-BODY SIMULATIONS 

In this section, we will briefly summarise the basic equa¬ 
tions that are used in f{R) cosmological simulations. 


A. Non-linear perturbation eqnations and screening 
mechanism 


as well as an equation for the scalar field fa = If the 
amplitude of the scalar field \fii \ is very small (|/rj| <C 1), the 
equation of motion for //j can be presented as 

^^SfR= ^[dR-SnGSp] , (7) 

where (j) denotes the gravitational potential, 5fa = fniR) — 
fniR), 5R = R—R, and dp = p—p. The overbar denotes the 
background quantities, and V is the derivative with respect to 


the physical coordinates. Combining Eq. (j^ 

and Eq. Q, we 

have 

= 47rG5p , 

(8) 

where 


A 

(PN — <P + 2 ’ 

(9) 


is the lensing potential ll25l . If we define an effective density 
field, the modified Poisson equation, Eq. (|^, can be cast into 
the same form as Eq. ([^ 

\7^^ = ATTGSpee , ( 10 ) 

where Jpeff = and G^g is the effective Newtonian 

constant which is defined by 

Geff characterizes the strength of gravity among massive par¬ 
ticles in f{R) gravity. 

As is well known, in linear theory f{R) gravity can be ruled 
out due to the enhancement of gravity relative to the stan¬ 
dard gravity Geff ~ |G. This conclusion is regardless of 
the functional form of f{R) (also see the appendix]^. Thus, 
a screening mechanism is essential for f{R) theories to evade 
stringent local tests of gravity. However, there are two aspects 
of the screening mechanism in f{R) theory to consider; 

• A high curvature, R ^ SirGp, should be recovered in 
high-density regions. 

• The fifth force, V//j, should be sufficiently suppressed 
in high-density regions as well. 

Given the fact that the functional form of f{R) is usually cho¬ 
sen as lim 7 j_>._|_oo I /fi I = 0, if a high curvature can be achieved 
in high-density regions, the fifth force will be automatically 
suppressed;V i5/h = V/^j —> 0. However, it is important 
to note that a high-density does not imply a high curvature 
in f{R) gravity. As is discussed in Ref. flEj . the condition 
for high curvature in high-density regions in f{R) gravity is 
closely related to the inequality 


The formation of large-scale structure in f{R) gravity is 
governed by the modified Poisson equation 
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Here, in addition to our previous numerical study ll251l . we also 
provide a strict mathematical proof of the above inequality. 
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making use of the Green’s function. The details can be found 
in Appendix [A| Based on Eq. ( [T2| i, a necessary condition for 
the high-curvature in high-density regions requires that the 
depth of the local potential well | — (/)| is close to or above the 
value of the background field | — (/)| > 2 c^|/ii;| Il25ll2^ : a suf¬ 
ficient condition for the low-curvature in high-density regions 
is that I — 01 <C 2c^\f]i\. If expressed in terms of the Newto¬ 
nian potential | — 0 Ar|, | — 0 Ar| > |c^|//j| is a weaker neces¬ 
sary condition for the high-curvature but | — <f>N\ <C |c^|/ij| 
is a stronger sufficient condition for the low-curvature in high- 
density regions. 


B. Cosmological simulations of f{R) gravity 

Our cosmological simulations are performed using the 
ECOSMOG code ll27l which itself is based on the A^-body 
code RAMSES 128]. The box size is Tbox = 150/i“^Mpc. 
The number of particles is iV = 256^. The cosmological pa¬ 
rameters are 0° = 0.049,0° = 0.267, = 0.684, h — 

0.671, ris = 0.962, and as = 0.834, which are the Planck 
2013 l3] best-fit values for the standard ACDM model. We 
use the Mpgrafic package OOl to generate initial conditions 
at z = 49. We choose the parameter fj^o for our f{R) model 
as /fio = -10 "®,/ro = -10"®,/ro = -10""^. In addition 
to the f{R) simulations, we also implement a suite of ACDM 
simulations with the same box size, cosmological parameters 
and initial conditions as control. 



“/fl 


FIG. 1. Upper panel: j3 = R{fR)/AK — 1 as a function of fu. It 
is clear that the accuracy of the fitting formula used in this work has 
been improved significantly particularly for small R. Lower panel: 
The error of fitting formulas. When R > 10i?o, the error is below 
1%. When R ~ Ro, the error is around 2%. The error is always less 
than 4% for R > 4(1 -|- l3m)R, where firn = 10"®. 


IV. POWER SPECTRUM 
A. The impact of background field fn 

As we have seen above, the scalar field /r enters the field 
equations Eq. (j^ and Eq. (j^ via i?(/R). Thus, to numerically 
solve the equations, we need to have analytical expressions 
of i?(/R), which can in principle be obtained by inverting 
fii{R), which itself can be derived from Eq. (j^. However, 
Eq.§ is complicated and it is difficult to extract R^fs) ana¬ 
lytically. Eollowing ESi we will instead use fitting formulae 

fori?(/R). 

In this work, we propose an improved fitting formula of 
Rifn) for our /(i?) model 


/3m = 10"®. 

To examine the impact of the accuracy of the fitting for 
RUr) on matter power spectra, we implement a test sim¬ 
ulation, using the same initial conditions as in our previous 
work Il26ll . The number of particles in the test is N = 256® 
and the box size is Lbox = 150/i"^Mpc. We choose /rq = 
— 10 "® for illustrative purposes since this model has the most 
complex shape of the power spectrum and most complicated 
screening situation at z = 0. The numerical results are shown 
in Eig. 1^ in which the upper panel shows the relative differ¬ 
ence between the power spectra of our /(i?) and the ACDM 
model 

AP/P = P^(ij)/PACDM - 1 , (14) 


P(/r) = 12^0772 + 


'Ir'' 


rje 


Hi , 


(13) 


where a and 77 are fitting parameters. Eig. shows the accu¬ 
racy of this improved fitting formula compared to the exact 
expression obtained from Eq. 0 and the one used in ||26|. 
p and a depend on 0 ^, and their values here are taken as 
rf = 1.3038, a — 0.3733 for (2^ = 0.316. However, 
they are independent of D and /rq. Eig. [T] shows that, when 
R > lOPo the error of the fitting formula is well below 1%, 
when R ^ Rq the error is around 2%, and the overall er¬ 
ror is always smaller than 4% for P > 4(1-1- /3m)A where 


at z = 0, measured by using the POWMES Il29l code. The 
lower panel shows the relative difference on the power spectra 
between this work and our previous work 


P 


previous work 


Pt 


this work 



X 100% 


We find a good agreement on the measured power spectra; 
the relative difference is below 1% up to fc ^ 5/i/Mpc, 
and even on the smallest scales probed by this simulation 
(k ^ 10/i/Mpc) it is less than 3%. We therefore conclude 
that the error induced by the fitting formula of P(/i^) has been 
controlled within a fairly reasonable range. 
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FIG. 2. The power spectra measured from our previous and current 
simulations using POWMES. The relative difference is well below 1% 
for k < 5/i/Mpc, and the difference on small scales k > 5/i/Mpc 
is less than 3%. 


B. Power spectra of scalar fields 


To achieve a high spatial resolution in simulations, the 
RAMSES code employs the adaptive mesh refinement (AMR) 
technique. The simulation box is initially covered by a do¬ 
main mesh with a fixed resolution. Each cell of the domain 
mesh is hierarchically refined according to some pre-defined 
criteria (e.g., a density threshold) during the simulation. Al¬ 
though the domain mesh resolution in our simulation is only 
256^, the highest resolution of the refined cells could be as 
high as (2^^)^ = 32768^. The physical quantities (e.g., den¬ 
sities, potentials) are assigned on the AMR grid, at the centers 
of the cells, during the simulation. 

In this work, we make use of this AMR grid to measure 
the power spectra of various scalar fields {5p, 5pes, 4 >n, 
c^5fii). Our method is similar to measuring the power spec¬ 
trum of gas pressure in hydrodynamical simulations Il28l . Un¬ 
like the density field, other scalar fields such as the gas pres¬ 
sure and various potentials can not be easily sampled by parti¬ 
cles without bias. The AMR grid is therefore a natural choice 
for this work. In order to do this, we record the values of these 
fields at different levels of the AMR grid for each snapshot. 
However, we only use the leaf cells (which are not refined) at 
each snapshot. The leaf cells can seamlessly cover the whole 
simulation domain. 

As discussed in Appendix [B] the power spectrum of a con¬ 
tinuous scalar field u{x), where u represents any one of 6p, 


6peff, 4>, (j>N and c^Sfa, is defined by 

(^UihWihr) = {27rf 6^{ki - k)Pu{k ), (15) 

where U{k) is the Fourier transform of u{x) and 5^ is the 
Dirac delta function. Although this definition assumes a con¬ 
tinuous scalar field, in practice the power spectrum can only 
be measured by discrete samplings. Usually, the scalar field 
is sampled on a set of regularly spaced grid points and then 
analyzed by using the Fast Fourier Transform (FFT). As ex¬ 
plained in Appendix the biggest problem of using FFT to 
measure the power spectrum is the aliasing effect, namely 
the discrete Fourier transform does not give the power spec¬ 
trum of the scalar field but a sum of aliases of the continuous 
Fourier transform of the scalar field 

= + , (16) 


where k^ = is the Nyquist frequency, L is the box size, 
Ng is the number of cells in one dimension of the FFT mesh, 
and n is a position vector which indicates the alias. Without 
prior knowledge of the power spectrum (e.g., its shape), it is 
impossible to remove the aliasing effect from the true signal. 
However, it is possible to minimize such aliasing effects by 
increasing the resolution of the FFT mesh. If the Nyquist fre¬ 
quency is high enough, the alias can be separated from the 
true signal in the main region 0 < k < 2kiy and so its effect 
can be minimized. 

In this work, we assign the values of scalar fields to the FFT 
grid without smoothing 


UXg - U{x)\g=Sg , 

where Xg represents the positions of the FFT grid points, and 
uSg is the sampled signal on those points. This is equivalent 
to interpolating the scalar field u{x) from the leaf AMR cells 
to the FFT grid points Xg directly. In detail: 

• If the leaf AMR cell is coarser than or the same coarse¬ 
ness as the FFT cell, we assign the value of the scalar 
field in the AMR cell to all of the FFT grid points that 
are contained within the AMR cell. 

• If the leaf AMR cell is finer than the FFT cell, we inter¬ 
polate values from the eight nearest surrounding AMR 
grid points to the FFT grid. 

This method enables us to measure the power spectra of any 
scalar field, and we call it the AMR-FFT method for conve¬ 
nience. Before presenting our results, we will perform several 
tests on this method. 

Our first test is about the accuracy of the AMR-FFT method 
when applied to the power spectrum of the dark matter den¬ 
sity field, Pjn{k). Pm{k) can be measured to a high accu¬ 
racy directly from the dark matter particles. We adopt the re¬ 
sults from the POWMES code as a control. POWMES is based 
on a method of Taylor expansion of trigonometric functions, 
and can yield an unbiased and nearly alias-free estimation of 
Pra{k). For the purpose of comparison, we only focus on one 
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realization, and choose f{R) models with different parame¬ 
ters as well as the ACDM model. These simulations share 
the same initial conditions. In Fig. we show the dark mat¬ 
ter power spectra for the ACDM model measured using our 
AMR-FFT method with different resolutions of the FFT grid. 
As shown in Fig. the accuracy of the AMR-FFT method 
depends strongly upon the resolution of the FFT grid. A low 
resolution measurement such as 512^ (fcAr/4 ~ 2.68/i/Mpc) 
gives very poor accuracy on the power spectrum on small 
scales k > kN/^- However, with a much higher resolu¬ 
tion, such as 4096^ (fcAr/4 ~ 20/i/Mpc), this method agrees 
with the POWMES code very well. Further, we also show 
in Fig. 0 the quantity AP/P as defined in Eq. (14i for dif¬ 
ferent f(R) models. Again, using the 4096^ FFT grid, the 
AMR-FFT method agrees with POWMES better than 1% out 
to fc ~ 10/i/Mpc. 


field as discrete mass-weighted dark matter particles 


Pm^H {x) 


^Pm{x) = ^J2mS°{x- Xi) 

i 

i 

'^mes{xi)6^{x - Xi) , 
i 


(18) 


in which m is the true mass of particles and Weff = is 

the effective mass. Although the discrete particles induce shot 
noise, the Fourier transform of the density field, in principle, 
can be evaluated accurately by a direct sum over the modes of 
the Fourier transform of each particle 


In Ref 133], a multi-grid method which is based on a hi¬ 
erarchy of nested cubic Cartesian grids is proposed in order 
to save the use of memory in the computation of the FFT. 
The power spectra can be measured by dividing the volume 
of a box into small cubes and then stacking the small cubes 
into a co-added density field. Instead of using a single high 
resolution FFT, only a few times of relatively low resolution 
FFT are needed in this method. The final power spectrum is 
obtained by combining the different “band power" measure¬ 
ments obtained from different volumes of the stacked density 
fields Il33l . In Fig. 0 we also present the power spectrum 
measured using this method. However, as shown in Fig.|^ the 
stacked density yields about 3% error on the power spectrum 
relative to POWMES. In order to get more accurate results, we 
therefore do not use this method in this work. 

Our second test concerns the measurement of the power 
spectrum of the effective density field, (5eff, which is defined 
by 


JT _ meff Geff 5p m 

Oeff = ^^- 


Using the AMR-FFT method, the power spectrum of <5eff can 
be measured in a similar way to that of the density field. For 
illustrative purposes, we only focus on the f{R) model with 
fjiQ — —10“®, and the results are shown in Fig. 0 As a 
comparison, we also present the dark matter power spectrum 
Pm{k). Compared with Pm{k), the effective power spec¬ 
trum Pm^ff{k) is enhanced due to the enhancement of grav¬ 
ity. However, Pm^ii{k) should not exceed ^Pm{k) since the 
maximal enhancement of gravity is 1 /3. This upper limit is 
indicated by the dashed line in Fig.0 

To further test the validity and accuracy of our AMR-FFT 
method on the power spectrum of the effective density field, 
instead of using the regular grid for sampling the effective 
density field, we sample the effective density field at the po¬ 
sitions of dark matter particles. We treat the effective density 


Pm^ajk) = y^^mes{xi)e . (19) 

i 

However, this method is numerically unfeasible for N-body 
simulations since they usually contain a large number of par¬ 
ticles and Eq. ([T^ does not have a fast algorithm like the EET 
due to the irregular distribution of points. In practice, in or¬ 
der to improve the efficiency, the power spectrum of discrete 
particles can be measured by assigning particles to a regular 
grid and then analyzing them using the EET. The assignment 
of particles is equivalent to smoothing the underlying density 
field and then putting the averaged values on an EET grid. This 
smoothing effect, however, can be exactly corrected after¬ 
wards. Furthermore, with the aid of further correction strate¬ 
gies (e.g. Ref Or] '), the aliasing effect can be significantly 
suppressed and a reasonable accuracy of the power spectrum 
can be obtained with less computational expense. Although 
this method is not fully free of bias and aliasing in general, 
it provides an alternative way to our AMR-FFT method for 
measuring the power spectrum of the effective density field. 
We therefore use it as a cross check of our measurement of 
p 

metf 

In high density regions 5^1, the mass weighted sam¬ 
pling should be unbiased and the power spectrum of the mass 
weighted particles should be close to the power spectrum 
of the effective density field. Since the POWMES code has 
not been tested for mass weighted particles, we instead use 
our own code. Our code follows the method as proposed in 
Ref 134|. We remove the aliasing effect by assuming a power 
law FWieff oc for the power spectrum on large k, where a 
is a fitting parameter. Following Ref. Il34ll . we also correct the 
shot noise induced by the discrete sampling of particles and 
the smoothing window function of the particle assignment. 
The results are shown as red stars in Fig. 0 As a demon¬ 
stration of our code, we also present the measurement of the 
dark matter power spectra Pm from the dark matter particles 
using our code. As is indicated by the blue stars, our code 
yields the same results as the POWMES code and the AMR- 
FFT method. As for the power spectrum of the effective den¬ 
sity field, on small scales the mass weighted sampling agrees 
very well with the results from our AMR-FFT method. How¬ 
ever, on large scales, there are some deviations. This can be 






6 



fc[/i/Mpc] 



FIG. 3. The power spectrum of the ACDM model measured from the 
POWMES code (the solid line) and our AMR-FFT method with dif¬ 
ferent resolutions. Our highest resolution measurement 4096^ agrees 
with the POWMES code very well out to fc ~ 10/i/Mpc. 


expected because the sampling of mass weighted dark mat¬ 
ter particles is biased in low density regions for the effective 
density field. In these regions, Eq. ( fTS] ! is less accurate and, 
moreover, there might be few or no dark matter particles in 
low density regions while the effective density field is contin¬ 
uous and can not be zero. 

In addition to the above tests, we also make several consis¬ 
tency tests. From the Poisson equations, Eq. ([^ and Eq. ( [T0] i, 
the power spectra of potentials and density fields should fol¬ 
low the relations 


FIG. 4. The relative differences of the conventional dark matter 
power spectra with respect to the ACDM model for different f{R) 
models. The 2048® FFT measurements yield about 5% relative dif¬ 
ferences on small scales (dotted lines), while the high resolution mea¬ 
surements 4096® (dashed lines) agree with POWMES (solid lines) 
better than 1%. The stars indicate the results obtained by stacking 
density fields. 


satisfy the following inequalities 

Pc^Sfnik) < ■^P4,{k) < . ( 21 ) 

InFig.|^ we show the ratio of \P^/Pc 2 sfj^ and ^P^j^/P^^gf 
as a function of k, from which we can see that indeed Eq. ( [2l| ) 
holds on scales from k ~ 0.06/i/Mpc to A: ~ 10/i/Mpc. 


P<pN 


- 0 ° 
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P^_p 

I J- rr. 




( 20 ) 


In our simulations, the potentials are obtained by using a re¬ 
laxation method in real space 03, which is different from the 
FFT method. It is therefore important to check the consistency 
of the above relations. In Fig. we show the directly mea¬ 
sured power spectra of </), (pj^ and the corresponding power 
spectra derived from the density fields (the right-hand sides of 
the above equations). Most strikingly, the numerical results 
agree very well over 12 orders of magnitude. This serves not 
only as a test of the consistency, but also as a check of the 
AMR-FFT method employed in this work. 

Finally, besides the above consistency relations, according 
to Eq. 0, the power spectra of the scalar fields should also 


C. Power spectra of the effective density field 

The above tests give us confidence in the validity and accu¬ 
racy of our AMR-FFT method. In this subsection, we present 
the power spectra of the effective density field measured by 
averaging over five realizations of simulations. Following the 
results of those tests, we use a 4096® FFT grid in this subsec¬ 
tion to measure the power spectra of density fields and various 
potentials. 

In Fig. we show the ratio of the power spectra of the two 
potentials 



Prr 


( 22 ) 


As indicated in Fig.l^ the ratio satisfies 1 < P^/P^^ < 16/9 
for all simulated /(it) models, consistent with the prediction 
by Eq. (21 1 . For fj^Q = —10“^, due to the lack of an efficient 
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FIG. 5. The power spectrum of the effective density field measured 
by our AMR-FFT method (red solid line). The red stars indicate the 
power spectrum of the mass weighted particles. For comparison, the 
power spectra of the dark matter density field are also shown. Our 
code (blue stars), the POWMES code (green triangles) and the AMR- 
FFT method (blue solid line) agree well on the power spectrum of the 
dark matter density field Pm- 
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FIG. 6. The power spectra of 0iv (blue triangles) and (p (red stars) 
measured using our AMR-FFT method as compared with the spectra 
derived from the density fields (solid lines). 
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FIG. 7. The ratio of \P^/Pc2sf,i and \P^^IP^i Sf^ as a function 
of k. 


Chameleon screening, the gravitational force is enhanced by 
1/3 and therefore is very close to 16/9. On the other 

hand, for //jq = —10“®, the screening mechanism works ef¬ 
ficiently and hence is very close to 1. The situation 

for fuQ = —10“^ is somewhere in between, as expected. 

Next, we show the relative difference of the effective power 
spectra with respect to the ACDM model. In Fig. we 
present AP„^/Pacdm and AP^^^j /Pacdm for the simulated 
f{R) models. As expected, compared with the dark mat¬ 
ter power spectra, the effective power spectra show more 
significant deviations from the ACDM model. For = 
— 10“^, APm^ff/PACDM is about three times as large as 
AP^/Pacdm. which is because the former contains contri¬ 
butions from both the latter and the fact that GeS is substan¬ 
tially enhanced compared with G. Even for fuQ = —10“®, 
the effective power spectra show sizeable deviations (about 
15%) from the ACDM model. 


D. The halo-halo power spectrum 

As shown in the previous subsection, the power spectra of 
P^j^ are quite different from P^ in f{R) gravity. This is an in¬ 
teresting feature, which can be used to discriminate the f{R) 
theory from other dark energy theories. To achieve this, how¬ 
ever, we need to measure the power spectra of P^^ and P^ 
independently. The lensing potential, can be measured di¬ 
rectly from weak lensing surveys. Upcoming galaxy surveys 
such as the Euclid mission iTSll have the power to measure 
P^^ accurately on large scales. 

The focus here is on how to measure P^, which is, indeed, 
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FIG. 8. The ratio of the power spectra = Pm^ff/Pm as a 

function of k. The lower and upper dashed lines indicate the value 
of 1 and 16/9, respectively. The ratio of the spectra should satisfy 
1 < P<t>/P<i>N < 16/9. For ~ —lO”"*, due to the lack of 
screening, Pt/j/Pr/,^ is very close to 16/9. For fno = —10~®, the 
screening mechanism works efficiently and so P^ /P<t>Ki is close to 1. 



fc[/i/Mpc] 


FIG. 9. The relative differences of power spectra with respect to the 
ACDM model for the dark matter density field (triangles with solid 
lines) and the effective density field (solid circles with dashed lines), 
respectively. The power spectra of the effective density fields 
show more significant deviations from the ACDM model than those 
of the dark matter density field Pm- 


nontrivial. One possible method is to measure the galaxy 
cluster-cluster power spectrum. In fact, from Fig. it can 
be seen that the differences between and also exist 
on relatively large scales. on large scales can be probed by 
the cluster-cluster power spectrum. In practice, it is more con¬ 
venient to work with the ratio Phh /Phh directly, where Phh 
is halo-halo power spectrum of standard halos. The advantage 
of using P™ /Phh is that the ratio is independent of the halo 
bias and is practically measurable. In observations, galaxy 
clusters can be observed using different methods such as x- 
ray observations, lensing and the thermal Sunyaev-Zeldovich 
effect 13^ . A wealth of information about the clusters can 
be obtained from these surveys. For example, the true mass 
can be inferred from the lensing data or the gas mass-cluster 
mass scaling relation El; the effective mass can be estimated 
either by the temperature-effective mass scaling relation lIMIl 
or by the profiles of gas density and temperature in x-ray sur¬ 
veys. In practice, we can first identify the clusters in x-ray 
surveys and measure the effective mass of each cluster. Then 
we divide the clusters into different mass bins. We measure 
Phh mass bin. Similarly, we can measure Phh for 

each mass bin as well. In the ACDM model, for a given mass 
bin the measurements of P™ and Phh should be the same 
Phh/Phh = 1- However, in /(P) gravity, Pf// and Phh can 
be different. There are two main reasons for the differences. 
The major one is due to mass calibration. Given a mass bin, 
the number of clusters with messes determined by effective 
mass will be quite different from the number determined by 
the true mass (see mass functions in Ref.lSl). The second 
reason is due to the positions of cluster centers. The centers 
of effective masses do not necessarily overlap with the true 
masses, especially for small clusters. Since the power spec¬ 
trum encompasses the above two effects, it should be more 
useful than the mass function to test f{R) gravity. 
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Next, we test the above idea using simulations. We con¬ 
struct halo catalogs using a modified version of the Amiga 
Halo Finder (Ahf) 091 . In halo-halo power spectrum Phh, 
halo number density fields are represented by un-weighted 
particles. We therefore can use the POWMES code to mea¬ 
sure the power spectrum. In Fig. we present the ratio of 
the effective halo power spectrum to the standard halo power 
spectrum for different f{R) models. Different colours repre¬ 
sent different mass bins. In Fig. 10 note that shot noise has al¬ 
ready been subtracted. On relatively large scales, the values of 
effective halo power spectra are less than those of the standard 
halos Phh /Phh < 1 for given mass bin. As shown in Fig. 
on relatively large scales {k < 0.2ft,/Mpc), for /^o = —10 
and f]io = —10“^, the ratio Pf/h /Phh deviates from the 
ACDM prediction (the dashed line) at a level of almost Icr. 
For ffio = —10“®, due to the screening mechanism, the ratio 
for massive clusters (M > is consistent with 

the ACDM prediction. However, for less massive clusters 
(M < Mq/ h), the deviations are over 3cr. Further, it 

can also be found that for //jo = —10“"^ and = —10“®, 
the errors on the ratio are greater than for Jbxi = —10“®. 
This is because the effective power spectra have larger scatters 
for the former two cases than for = —10“® as shown in 
Fig. El The large errors on small scales (fc > 0.2/i/Mpc) are 
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fm^ 10 ^ fm^ 10 ® fm^ 10 



k[h/Mpc] k[h/Mpc] k[h/Mpc] 


FIG. 10. The ratio of the effective halo power spectrum to the standard halo power spectrum for different f{R) models. is measured 
from the effective catalog and P^h is measured from the standard catalog. Different colours represent different mass bins. In effective catalog, 
effective mass of halo is used and in standard halo catalog, the conventional mass is used. On relatively large scales (fc < 0.2/i/Mpc), for 
fuo = —10“'^ and /ijo = —10“®, the ratio deviates from the ACDM prediction (the dashed line) at a level of Icr. For /jjo = —10“®, due 
to the screening mechanism, the ratio for massive clusters (M > is consistent with the ACDM prediction. However, for less 

massive clusters (M < the deviations are over Scr. 


due to shot noise. In particular, for the most massive clusters 
(M > we only have a few hundred samples. 

The level of shot noise is quite high there. Finally, due to the 
limited box size and number of realizations, the simulations 
performed in this work tend to over estimate the errors. With 
a larger simulation box and more realizations, the error bars 
can be reduced and the deviations of f{R) models from the 
ACDM model should be more clear. The results presented 
above therefore are very conservative. 


V. SUMMARY AND DISCUSSION 


In this work, we have studied the non-linear power spectra 
of scalar fields in f{R) gravity, using a suite of N-body sim¬ 
ulations. We have illustrated in detail our AMR-FFT method 
for measuring the power spectra of scalar fields. Using this 
method, we have measured the power spectra of the potentials 
and the effective density fields. Our main results are summa¬ 
rized as follows: 


• We have verified the inequality 


c" \5fR\ < 




(23) 


in Fourier space by comparing the power spectra of the 
potentials. As is discussed in Ref ll26l . the above in¬ 
equality is closely related to the screening mechanism 
in f{R) gravity. Its validity is important for predicting 
the screening theoretically. 


• We find that compared with the dark matter power spec¬ 
tra, the effective power spectra differ more significantly 
from the ACDM model. Even for fuQ = —10“®, the 
effective power spectrum shows a sizeable signal of de¬ 
viation. 

We have tested that these conclusions are applicable to other 
f{R) models as well. For Hu-Sawicki model ESi . we find 
similar results as presented in Ref. iQ]. 

Since the formation and clustering of galaxies in f{R) 
gravity is more closely related to the effective density field 
rather than the dark matter density field itself, our findings in¬ 
dicate that the traditional statistics of the dark matter density 
field such as the power spectrum or equivalently the two-point 
correlation function may seriously under-predict the impact of 
modifications of gravity on the clustering of galaxies. How¬ 
ever, it should be noted that although the power spectrum of 
the effective density field in f{R) gravity is significantly dif¬ 
ferent from that in the ACDM model, we should caution that 
these large deviations are likely to be highly degenerate with 
galaxy bias. Robust constraints on f{R) gravity can thus only 
be obtained once we have a solid knowledge of the expected 
biasing of galaxies in f{R) cosmologies. 

Nevertheless, it is interesting to find that the power spec¬ 
trum is quite different from in f{R) gravity. This can 
be used to discriminate the theory from the ACDM model. 
We tested this idea by investigating the halo-halo power spec¬ 
trum on relatively large scales. Galaxy clusters can be ob¬ 
served using different methods and a wealth of information 
about the clusters can be obtained. We found that probing 
the ratio P^h /Phh is an useful way to test f{R) gravity. On 
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relatively large scales (fc < 0.2/i/Mpc), for fjiQ = —10“"^ 
and fjio = —10“^, the ratio deviates from the ACDM pre¬ 
diction at a level of almost Icr. For /jjq = —10“®, due 
to the screening mechanism, the ratio for massive clusters 
(M > IO^^-^Mq/Zi) is consistent with the ACDM prediction. 
However, for less massive clusters {M < Mq/ h), the 

deviations are over Scr. In fact, due to the limited box size and 
number of realizations, the simulations performed in this work 
tend to over estimate the errors. Our results therefore can be 
further improved with a larger simulation box and more real¬ 
izations. The results quoted the above are very conservative. 
It is interesting to note that the latest all-sky Planck catalog 
of Sunyaev-Zeldovich (SZ) sources BD has already accu¬ 
mulated over 10^ clusters. The clusters are distributed over a 
large area of the sky. If the pointed X-ray follow-up informa¬ 
tion is available in the future, the catalog can be used to test 
f{R) gravity. Further, upcoming the eROSITA survey 1421 
will have the ability to discover over 10^ clusters. The shot 
noise of the measured power spectrum can be reduced signif¬ 
icantly. 

Our idea can also be further extended. It is interesting to 
know that the ratio is strongly scale dependent es¬ 

pecially when k < 0.2Zi/Mpc (see. Fig. [^, on which scales 
the galaxy bias is expected to be approximately scale indepen¬ 
dent. Unlike the halo number density field, the total number 
density field of galaxies, including both central and satellite 
galaxies, ought to have a close tie with the effective density 
field. This is because, for massive halos M > 
the total number of satellite galaxies formed inside such halos 
should relate to the effective mass rather than the true mass of 
the halos. Thus, if we work with the assumption that galax¬ 
ies are tracers of the effective density field then their clus¬ 
tering properties can be used to infer the power spectrum of 
the potential P^. The power spectrum of the lensing poten¬ 
tial can then be inferred from weak lensing data. Future 
galaxy surveys such as the Euclid mission llTSll will measure 
both of these statistics in the same region of the sky on large 
scales. By combining these two pieces of information we may 
be able to place robust constraints on f{R) gravity. However, 
this method requires an understanding of galaxy formation in 
f{R) gravity. We therefore shall elaborate on this idea in de¬ 
tail in future work. 
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Appendix A: Inequalities 


it to scalar fields. 

In the linear regime, 5R can be linearised as 

5R « -^5fR = 3c^mlg6fR , (Al) 

JRR 


where Jrr) is the effective mass for the scalar 

field. Eor a point mass particle, under the boundary conditions 
limr_>.+oo = lim,—i.+oo = 0, the linearized Eq.Q 

and Eq.([^ have solutions 


SfR 


2Gm 
3c^ r 
Gm 


(A2) 


where m is the mass of the particle. The potential SJr is of 
the Yukawa type which will be suppressed on scales above the 
Compton length Ac = l/rrieff. Erom Eq.lj^, we obtain 


Gm 

r 


Gm 


(A3) 


On large scales r ^ Ac, therefore. 


g-meffr ^ Q 


the gravity (p goes back to the standard gravity ~ (^jv- How¬ 
ever, on small scales r <C Ac 


g-meffr ^ ^ ^ 

f{R) gravity has a 1/3 enhancement relative to standard grav¬ 
ity (j) ^ /(P) gravity is therefore invalid in the linear 

regime regardless of the functional form of /(P). 

Eurther, in the large-held limit 131113^ . meff —0, the ex¬ 
ponential term in Eq.( |A2| l becomes 

g-meffr ^ ^ 


obtains its maximal value, 
treme case, we have 


Erom Eq.( A3 i, in this ex- 


C^\SfR\ 



(A4) 


On the other hand, in the small-held limit OTIl . mes +oo 
(e.g. the early Universe, when perturbations are small) , the 
potential S/r will be signihcantly suppressed 


ISJrI < 



(A5) 


In general cases, for the linearized equations, given a hnite 
volume V of the density held, under the boundary condition 

lim u{x) = 0 , 

|;e|—>-+ oo 


the potential is simply the superposition of the potentials gen¬ 
erated by local density helds. 


In this appendix, we provide a strict proof of Eq. ( [T^ . We 
shall start with the case of a single particle and then generalize 


u{x) 


dx' G{x, x')6p{x') 


(A6) 


V 
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where u{x) standards for the scalar fields and (j), 

respectively. G{x,x') is the Green’s function which is given 
by 


G{x,x') = <- 


2G 

Z\x—x' I 

G 


\x—x' I 

G 

\x—x' I 


G ^-m^.iAx-x\ 
3|^— 


icHfi,), 

(4>n), 

(0). 

(A7) 


In high-density regions Sp ^ 1, the contribution from the 
low-density regions {Sp < 0) to the total scalar field compared 
to the local contribution from the high-density region itself, 
can be neglected. From the Green’s functions Eq. ( |A7[ ), it 
follows that 


c" \5fR\ < 




(A8) 


The above inequality holds everywhere in high-density re¬ 
gions |l26l. However, in low-density regions, the inequality 
Eq.([T2)i may not hold everywhere because the local contribu¬ 
tion from the under-dense regions {6p < 0) may cancel out 
the contribution from distant high-density regions to the local 
total scalar field. Sfn and —(j)N may have different signs. The 
absolute values of 5Jr and —(pN may have very complicated 
relations. 


Appendix B: Alias sum 


Let u{x) be a continuous scalar field. The Eourier transform 
and its inverse transform are given by 


U(k) = j (fxu{x) 
(fik 


i{x) = 


(27r)' 


C/(fc) 


ik-: 


The two-point correlation function is defined by 

{u{xi)u{x2)*) 


dPkid^k2 

(2^)6 


U{ki)U{k2T)e 


*\ pi{kl-Si-k2-X2) 


(Bl) 


(B2) 


If we define the power spectrum of the scalar field u{x) as 
' t/(fcl)C/(fc2)*) = {2TTfS{ki - k2)Pu{k) , (B3) 

then the two point correlation function can be written as 


with n^jUy, Uz being integers. On the discrete grids, the in¬ 
tegration of Eq. •El) can be approximated as a sum over cells 
with volume dx^ « ^ 


S n' 


-ik-Xn 


(B5) 


where = u{x)\g=g^. The inverse Eourier transform of 
Eq. ( |B5[ ) can be presented as 

l■2kJs, ^2fcjv /■2fcjv Jj, Rh Rl. _ - 

= I I I , (B6) 


/o ^0 ^0 


(2 


TT 


where k^ = is the Nyquist frequency. The relation be¬ 
tween U(k) and U(k) can be obtained by noting that 






U{k)P^-^^ 


^ (27r)3 

f.2feiv /■2fcjv /■2fciv jj. Ri, jjL __ 


/O ^0 ^0 


where in the last equality we have used 


(B7) 


^ik-Xg _ ^i{k-\-2kNh)-: 


(B8) 


Comparing Eq. (B7 1 with Eq. 


we obtain 


u{k) = U{k + 2kNn) 


(B9) 


The discrete Eourier transform is simply the sum of replicas 
of the continuous Fourier transform. This result is known as 
the “alias sum”. 


Next, we evaluate the two point correlation Eq. ( ^ i on dis¬ 
crete grids. Given that dkf = dfcf « (x)^’ ^9- ' 
approximated as 


{u{xi)u{X2)*) 

ii E 


k2 


^i(ki-xi-k2-X2) 


ki,k2 




(27r)3 ^ V T 

k 


X 2 ) 


L3 


, ^ f j;,3 (Efcl } ik-(xi-X2) 

'(2^)3 J 7.3^ 


(BIO) 




X 2 ) = {u{xi)u{x2)*) 


1 

(27r)3 




where k = ‘^n indicates the discrete grids in the Fourier 
(B4) space and = U{k)\^^f,. In the above derivations, we have 
used 


In practice, we can only treat the continuous scalar field 
u{x) on discrete grids Xg = , where Ng is the num¬ 

ber of grid cells in one dimension, L is the box size, and 
n = UxX -f nyjj -f UzZ describes the position of grid points 


so that the correlation function is only dependent on the spa- 
cial separation r = £2 — xi. 
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Comparing Eq. (BIO i with Eq. ( |B4[ ), we have 


Pu{k) = 


L3 


Using Eq. (B91, we have 


where FFT[rtjg] = is the fast Fourier trans¬ 

form of the discrete grid points . 

Under the assumption of ergodicity, the ensemble average 
can be replaced by a spatial average. The isotropic power 
spectrum can be estimated by 


® (B13) 

Sn (|f^fc+2fejvnl^) 

^ L3 —- 

= '^Puik+ 2kNn) , 
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